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Abstract 



Heat conduction of one-dimensional chain of equivalent rigid particles in the field of external on- 
site potential is considered. Zero diameters of the particles correspond to exactly integrable case with 
*73 ' divergent heat conduction coefficient. By means of simple analytical model it is demonstrated that for 

any nonzero particle size the integrability is violated and the heat conduction coefficient converges. The 
q | result of the analytical computation is verified by means of numerical simulation in a plausible diapason 

of parameters and good agreement is observed 

c3 1 Introduction 

Heat conductivity in one-dimensional (ID) lattices is a well known classical problem related to the micro- 
scopic foundation of Fourier's law. The problem started from the famous work of Fermi, Pasta, and Ulam 
(FPU) [1], where an abnormal process of heat transfer was detected for the first time. Non-integrability of a 
system is a necessary condition for normal heat conductivity. As it was demonstrated recently for the FPU 
lattice [2, 3, 4], disordered harmonic chain [5, 6, 7], diatomic ID gas of colliding particles [8, 9, 10, 11], and 
the diatomic Toda lattice [12], non-integrability is not sufficient in order to get normal heat conductivity. It 
leads to linear distribution of temperature along the chain for small gradient, but the value of heat flux is 
proportional to 1 /N a , where N is the number of particles of the chain, and number exponent < a < 1. 
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Thus, the coefficient of heat conductivity diverges in the thermodynamic limit N — > oo. Analytical estima- 
tions [4] have demonstrated that any chain possessing an acoustic phonon branch should have infinite heat 
conductivity in the limit of low temperatures. 

From the other side, there are some artificial systems with on-site potential (ding-a-ling and related 
models) which have normal heat conductivity [13, 14]. The heat conductivity of the Frenkel-Kontorova 
chain was first considered in Ref. [15]. Finite heat conductivity for certain parameters was observed for 
Frenkel-Kontorova chain [16], for the chain with sinh-Gordon on-site potential [17], and for the chain with 
4 on-site potential [18, 19]. These models are not invariant with respect to translation and the momentum is 
not conserved. It was supposed that the on-site potential is extremely significant for normal heat conduction 
[18] and that the anharmonicity of the on-site potential is sufficient to ensure the validity of Fourier's law 
[20]. A recent detailed review of the problem is presented in Ref. [21]. 

Probably the most interesting question related to heat conductivity of ID models (which actually in- 
spired the first investigation of Fermi, Pasta and Ulam [1]) is whether small perturbation of an integrable 
model will lead to convergent heat conduction coefficient. One supposes that for the one-dimensional chains 
with conserved momentum the answer is negative [22]. Still, normal heat conduction has been observed in 
some special systems with conserved momentum [23, 24, 25], but it may be clearly demonstrated only well 
apart from integrable limit. 

The situation is not so clear in the systems with on-site potential. Although it was supposed that non- 
integrable system without additional integrals of motion would have convergent heat conductivity [22], no 
rigorous proof was presented. From the other side, recent attempt of numerical simulation of heat transfer 
in Frenkel-Kontorova model [26] demonstrated that because of computational difficulties no unambiguous 
conclusion can be drawn whether heat conduction is convergent for all finite values of perturbation of the 
integrable limit system (linear chain or continuous sine-Gordon system). 

It seems that computational difficulties of investigation of heat conduction in a vicinity of integrable 
limit are not just issue of weak computers or ineffective procedures. In the systems with conserved momen- 
tum divergent heat conduction is fixed by power-like decrease of heat flux autocorrelation function with 
power less than unity. Still, for the systems with on-site potential exponential decrease is more typical [26]. 
For any fixed value of the exponent the heat conduction converges; if the exponent tends to zero with the 
value of the perturbation of the integrable case, then for any finite value of the perturbation the character- 
istic correlation time and length will be finite but may become very large. Consequently, they will exceed 
any available computation time or size of the system and still no conclusion on the convergence of heat 
conduction will be possible. 

In this respect it is reasonable to mention the findings of our papers [24, 25, 26] concerning the transi- 
tions from infinite to finite heat conduction in the chain of rotators and Frenkel - Kontorova system. These 
findings were criticized in paper [27] by providing the computation results for larger systems. We agree 
with the authors that the results of papers [24, 25] do not prove the reality of the above transition. In fact, we 
have pointed it out in [26] . From the other side, it is clear that such more extended numerical simulations 
cannot prove that the transitions do not exist at all. Probably, such conclusion may be driven if, with the 
help of some not yet developed method, the simulation length or effective time will be extended by many 
orders of magnitude. 

The other way to overcome this difficulty is to construct a model, which will be, at least to some extent, 
tractable analytically and will allow one to predict some characteristic features of the heat transfer process 
and the behavior of the heat conduction coefficient. Afterwards the numerical simulation may be used to 
verify the assumptions made in the analytic treatment. To the best of our knowledge, no models besides 
pure harmonic chains were treated in such a way to date. Introduction and investigation of a nontrivial 
model of this sort is a scope of present paper. 

We are going to demonstrate that there exist models which have integrable system as their natural 
limit case, small perturbation of the integrability immediately leading to convergent heat conduction. The 
mechanism of energy scattering in this kind of systems is universal for any temperature and set of the 
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Figure 1 : Sketch of the hard-disk chain exposed to periodic on-site potential U (x) (a is the period of the potential, 
Uq - its height, d is the diameter of the disks). Piecewise-linear potential (10) (a) and sinusoidal potential (20) (b) are 
plotted. 

model parameters. The simplest example of such model is one-dimensional set of equal rigid particles with 
nonzero diameter (d > 0) subjected to periodic on-site potential. This system is completely integrable only 
if d = 0. It will be demonstrated that any d > leads to effective mixing due to unequal exchange of energy 
between the particles in each collision. This mixing leads to diffusive mechanism of the heat transport and, 
subsequently, to convergent heat conduction. 



2 Description of the model 

Let us consider the one-dimensional system of hard particles with equal masses subject to periodic on-site 
potential. The Hamiltonian of this system will read 

H = Y.{\ M *l + V(Xn + l ~ X n ) + U( Xn )}, (1) 

where M - mass of the particle, x n - coordinate of the center of the n-th particle, x n - velocity of this 
particle, U (x) - periodic on-site potential with period a [U (x) = U(x + a)]. Interaction of absolutely hard 
particles is described by the following potential 

V(r) = oo if r < d and V(r) =0ifr > d, (2) 

where d is the diameter of the particle. This potential corresponds to pure elastic impact with unit recovery 
coefficient. Sketch of the model considered is presented at Fig. 1. 

It is well-known that the elastic collision of two equal particles with collinear velocity vectors leads to 
exchange of their velocities. An external potential present does not change this fact, since the collision takes 
zero time and thus the effect of the external force on the energy and momentum conservation is absent. 

The one-dimensional chain of equivalent hard particles without external potential is a paradigm of the 
integrable nonlinear chain, since all interactions are reduced to exchange of velocities. In other words, 
the individual values of velocities are preserved and just transferred from particle to particle. It is natural 
therefore to introduce quasiparticles associated with these individual values of velocities. They will be char- 
acterized by a pair of parameters (E k , n fc ), where E k = v\j2 is an energy of the quasiparticle, n fc is a unit 
vector in a direction of its motion. Every particle in every moment "carries" one quasiparticle. The elas- 
tic collision between the particles leads to simple exchange of parameters of the associated quasiparticles, 
therefore the quasiparticles themselves should be considered as non-interacting. 
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3 Analytical study 



The situation changes if the external on-site potential is present. It is easy to introduce similar quasi- 
particles (Ek will be a sum of kinetic and potential energy). The unit vector n of each quasiparticle between 
subsequent interactions may be either constant (motion in one direction) or periodically changing (vibration 
of the particle in a potential well). In every collision the particles exchange their velocity vectors, but do 
not change their positions. Consequently two quasiparticles interact in a way described by the following 
relationships: 

E[ = E 1 + U(x c + d/2)-U(x c -d/2), 
E' 2 = E 2 -U(x c + d/2) + U(x c -d/2), 

ni = m, (3) 
n' 2 = n 2 . 

The values denoted by apostrophe correspond to the state after the collision, x c is a point of contact between 
the particles. It should be mentioned that in the case of nonzero diameter the quasiparticles are associated 
with the centers of the carrying particles. 

If the diameter of the particles is zero, then the additives to the energies in the first two equations of 
system (3) compensate each other and the energies of the quasiparticles are preserved in the collision. 
Therefore effectively the interaction between the quasiparticles disappears and the chain of equal particles 
with zero size subject to any on-site potential turns out to be completely integrable system. Thus, contrary 
to some previous statements it is possible to construct an example of strongly nonlinear one-dimensional 
chain without momentum conservation, which will have clearly divergent heat conductivity (even linear 
temperature profile will not be formed). 

The situation differs if the size of the particles is not zero, as the individual energies of the quasiparticles 
are not preserved in the collisions. In order to consider the effect of such interaction we propose simplified 
semi-phenomenological analytical model. 

After I collisions the energy of the quasiparticle will be 

i 

E(l) = E + AEj, AEj = U(xj + d/2) - U{xj - d/2), (4) 
j'=i 

where j-th collision takes place in point Xj, E is the initial energy of the quasiparticle. Now we suppose 
that the coordinates of subsequent contact points {..., Xj,Xj + i, ...}, taken by modulus of the period of 
the on-site potential, are not correlated. Such proposition is equivalent to fast phase mixing in a system close 
to integrable one and is well-known in various kinetic problems [28]. The consequences of this proposition 
will be verified below by direct numerical simulation. 

Average energy of the quasiparticle is equal to (E ) over the ensemble of the quasiparticles, as obviously 
(AEj) = 0. Still, the second momentum will be nonzero: 

((£(/) - E ) 2 ) = l{(U(x + d/2) - U(x - d/2)) 2 ) x (5) 

The right-hand side of this expression will depend only on the exact shape of the potential function 

((E(l) - E ) 2 ) = lF(d), F(d) = - f a [U(x + d/2) - U(x - d/2)] 2 dx, (6) 

a Jo 

The last expression is correct only at the limit of high temperatures; it neglects the fact that the quasiparticle 
spends more time near the top of the potential barrier due to lower velocity. 

Let us consider the quasiparticle with initial energy E > Uq, where Uq is the height of the potential 
barrier. Therefore vector n is constant. Equations (5) and (6) describe random walks of the energy of the 
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quasiparticle along the energy scale axis. Therefore after certain number of steps (collisions) the energy 
of the quasiparticle will enter the zone below the potential barrier E(l) < U . In this case the behaviour 
of the quasiparticle will change, as constant vector n will become oscillating, as described above. After 
some additional collisions the energy will again exceed Uq, but the direction of motion of the quasiparticle 
will be arbitrary. It means that the only mechanism of energy transfer in the system under consideration 
is associated with diffusion of the quasiparticles, which are trapped by the on-site potential and afterwards 
released in arbitrary direction. Such traps-and-releases resemble Umklapp processes of phonon-phonon 
interaction [28], but occur in a purely classic system. 

The diffusion of the quasiparticles in the chain is characterized by mean free path, which may be evalu- 
ated as 

. 2a((U - E ) 2 ) 2a[2(k B T) 2 -2U k B T + U$] 

n c F(d) n c F(d) ' V ; 

where n c is a number of particles over one period of the on-site potential (concentration). Coefficient 2 
appears due to equivalent probability of positive and negative energy shift in any collision, T - temperature 
of the system, k B - Boltzmann constant. 

Average absolute velocity of the quasiparticle may be estimated as 



,. .. a irkfiT 

M A— — (8) 

Here the first multiplier takes into account the nonzero value of d and absolute rigidity of the particles. The 
second one is due to standard Maxwell distribution function for ID case. 

Heat capacity of the system over one particle is unity, as the number of degrees of freedom (i.e. the 
number of quasiparticles) coincides with the number of the particles and does not depend on the temperature 
and other parameters of the system. Therefore the coefficient of heat conductivity may be estimated [28] as 

2a 2 2(k B T) 2 - 2U k B T + Uf 
n c (a — n c d) F(d) 



«~A(H) zi^l "T"^ ' ° Jirk B T/2 (9) 



It is already possible to conclude that according to (9) regardless the concrete shape of potential U(x) 
in the limit d — > we have F(d) — > and therefore n — > oo, although for every nonzero value d the 
heat conductivity will be finite. Therefore unlikely known models with conserved momentum the small 
perturbation of the integrable case d — immediately brings about convergent heat conductivity. 



4 Numerical simulations 

It is convenient to introduce the dimensionless variables for the following numerical simulation. Let us set 
the mass of each particle M — 1, on-site potential period a = 2, its height U = 1, and Boltzmann constant 
k B = 1 in all above relationships. We suppose that the chain contains one particle per each period of the 
potential, i.e. that n c = 1, and the particle diameter < d < 2. 
Let us consider periodic piecewise linear on-site potential 

U(x) = x if x e [0, 1], 

U(x) = 2-x if x E [1,2], (10) 
U(x + 2l) = U(x) for x 6 [0,2], I = 0, ±1, ±2, ... 

(the shape of the potential is presented at Fig. 1). Then it follows from (9) that the non-dimensional heat 
conduction coefficient is expressed as 

8(2T 2 - 2T + 1) 



(2 - d)F(d) 



vrT/2, (11) 
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where function 



F(d) = d 2 -2/3d 3 , forO<d<l (12) 
F(d) = -4/3 + 4d-3d 2 + 2/3d 3 , for 1 < d < 2. 

The numerical scheme for solving the equations of motion describing the dynamics of the ID hard-point 
gas has been developed in a series of papers [29, 8, 30]. 

Dynamics of the system of particles with potential of the nearest-neighbor interaction (2) and piece- 
wise linear on-site potential (10) may be described exactly. Between the collisions each particle moves 
under constant force with sign dependent on the position of the particle. Therefore the coordinate of each 
particle depends on time t as piecewise parabolic function which may be easily computed analytically. If the 
particle centers are situated at distance equal to d, then elastic collision occurs. The particles exchange their 
momenta as described above and afterwards the particle motion is again described by piecewise parabolic 
functions until the next collision. 

Let us consider finite chain of N particles with periodic boundary conditions. Let at the moment t = 
one particle be at each potential minimum and let us choose Boltzmann's distribution of the initial velocity. 
Solving the equations of motion, we find a time t\ of the first collision between some pair of the adjacent 
particles, next a time t 2 of the second collision, in general between another pair of the adjacent particles, and 
so on. As a result, we obtain a sequence {U, where U is the time of the zth collision in the system, 

and Ui and rij + 1 are the particles participating in this collision. Since we need to implement numerical 
simulations as long as possible, in order to find the time asymptotic of the heat flux autocorrelation function 
entering the Green-Kubo formula, we use the numerical scheme of paper [9]. First, we incorporate the 
energy change of the n^th particle during the ith collision as 

= \{ V 'th ~ V *h) = \{Vn l+ l 2 ~ V ni 2 ). 

Next, we introduce a time step At, which is significantly less than the simulation time, but satisfies the 
inequality At t , where to = lim i ^ 00 (tj/z) is the mean time between successive collisions. Then, for 
each k = 0, 1, we define the local energy flow as a piecewise constant (in time) function 

3n{t) = ^ E AE n„ kAt<t< (k + l)At, (13) 

where the integer sets ifc n 's are defined by 

I kn = {i\ kAt < U <(k + l)A t , rii = n}. 

The set I kn takes into account those collisions that occur between particles n and n + 1 during the time 
interval kAt < t < (k + 1) At. Equilibration times were typically occurring in the system of the order 10 6 . 
After these times have passed, we define the time-averaged local thermal flow 

J n = (j n (t)) t = lim - f j n (r)dr (14) 

t->co t JO 

and the temperature distribution T n = (v 2 (t))t, where v n (i) is the velocity of particle n calculated at a time 
t. To find these averaged quantities, we have used times up to 10 7 . 

To find the heat flux autocorrelation function C(t) numerically, we calculated the time mean 
(J(t)J(t — t)) T /NT 2 , with J(t) = ^2 n j n (t) being the total heat flow through the gas/chain system con- 
sisting of N = 500 particles and temperature T = X) n T n /iV averaged over 10 4 realizations of initial 
thermalization. 

Numerical simulation of the dynamics demonstrates exponential decrease of the autocorrelation C (t) ~ 
exp(— at) for all values of the diameter < d < 2 and temperature T > where the simulation time is 
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Figure 2: CoiTelation function of the system of particles with d 
(curves 1, 2 and 3). 



0.5 under temperatures T = 0.24, 0.45 and 0.75 



plausible from technical viewpoint. For low temperatures however the exponential decrease is accompanied 
by oscillations with period corresponding to the frequency of the vibrations near the potential minima (Fig. 
2). The reason is that if the temperatures are low, the concentration of transient particles decreases expo- 
nentially and majority of the particles vibrates near the potential minima. It means that the ID gas on the 
on-site potential has finite heat conductivity. Coefficient of the exponential decrease of the autocorrelation 
function 

lnC(t) 



a 



lim 

t— >oo 



t 



and coefficient of the heat conduction 



K 



C(t)dt. 



(15) 



(16) 



are computed numerically. 

Dependence of a and k on particle diameter d is presented at Fig. 3. Maximum of a and minimum of k 
is attained at d = 1.4. As the temperature grows, a decreases [Fig. 4 (a)], and heat conduction k increases 
[Fig. 4 (a)]. 

Theoretical analysis of the heat conductivity presented above allows only approximate [although rather 
reliable, see Figs. 3(b), 5] prediction of the numerical value of the heat conduction coefficient k. Still, 
the other question of interest is the asymptotic dependence of the heat conduction on the parameters of the 
model. Formulae (11), (12) lead to the following estimations: 



K ~ 



K 



rp5/2 

<d-\ 



for T 
for d - 



oo, 
+0, 



(17) 
(18) 
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Figure 3: Dependence of the coefficient of the exponential decrease of the autocorrelation function a (a) and the 
coefficient of the heat conduction k (b) on the particle diameter d of ID gas at T=l. Curves 1 and 2 correspond to 
piecewise linear on-site potential (10), curve 3 represent theoretical predictions according to formula (11). 
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Figure 4: Temperature dependence of exponent coefficient a (a) and heat conduction coefficient k (b) for particle 
diameter d = 0.5. 
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Figure 5: Dependence of heat conduction coefficient on the temperature. The markers correspond to numerical 
results (In k versus In T), the straight line is In k = 2.5 In T + 3.45, corresponding to ( 17). Particle diameter d = 0.5 . 

k ~ (2 - cfT 3 , for d -> 2 - 0. (19) 

These estimations should be compared to numerical results. 

In order to check estimation (17) we consider the dependence of the logarithm of the heat conduction 
In k on the logarithm of the temperature In T. From Fig. 5 it is clear that in accordance with (17) In k grows 
as 2.5 In T as T — > oo. Fig. 6 (a) demonstrates that as d — > +0, the logarithm InK grows as — 2 hid, in 
accordance with (18). Fig. 6 (b) demonstrates that as d — * 2 — 0, the logarithm In k grows as —3 ln(2 — d), 
in accordance with (19). So, it is possible to conclude that analytical estimations (17), (18) and (19) fairly 
correspond to the numerical simulations data. 

The above analytical estimations imply that the type of dependence of characteristic exponent a and 
heat conductivity k on diameter d and temperature T does not depend on the concrete shape of on-site 
potential U (x) - actually, only its finiteness and periodicity do matter. Piecewise linear periodic potential 
(10) was chosen since it allowed essential simplification of the numerical procedure. For comparison we 
have considered also the smooth sinusoidal periodic potential 

U(x) = [1 -cos(ttw)]/2 (20) 

with period 2 and amplitude Uq = 1, similarly to potential (10). 
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In d 




-1.5 -1 -0.5 0.5 



In (d-2) 

Figure 6: Dependence of the heat conduction coefficient on the particle diameter (logarithmic coordinates, In k versus 
In d (a) and versus ln(2 — d) (b), curves 1 and 3). Lines In k = — 2 In d (curve 2) and In k = —3 ln(2 — d) (curve 4) 
correspond to relationships (18) and (19). Temperature T = 1 . 
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Figure 7: Dependence of the coefficient of the exponential decrease of the autocorrelation function a (a) and the 
coefficient of the heat conduction k (b) on the particle diameter d of ID gas at T=l. Curves 1 and 2 correspond to 
smooth on-site potential (20), curve 3 represent theoretical predictions according to formula (11). 
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Potential (20) does not allow exact integration and requires standard numerical procedures. Therefore 
it is convenient also to change rigid wall potential (2) by smooth Lennard- Jones potential 



Parameter e > characterizes rigidity of the potential, the hard-particle potential being the limit case: 



Methods of computing of the autocorrelation function C{t) and the heat conduction coefficient k in ID 
chain with analytic potentials of interaction are described in [26] in detail. It should be mentioned that in 
order to get close to the limit of the hard particles we should use small values of e (at the temperature T = 1 
value e = 0.01 was used). It implies rather small value of the integration step. (We used standard Runge- 
Kutta procedure of the fourth order with constant integration step At = 0.0001). Therefore for the case of 
hard (or nearly hard) particles the simulation with smooth on-site potential (20) is far more time-consuming 
than the simulation with piecewise linear potential (2). 

In the case of hard particles with smooth on-site potential the autocorrelation function C(t) decreases 
exponentially as t — > oo for all range < d < 2, T > 0, i.e. the heat conduction converges. Fig. 
7 demonstrates that the type of dependence of a and k on parameters d and T is similar for piecewise 
linear potential (2) and sinusoidal potential (20) (although numerical values a and k vary slightly). For 
this potential function F(d) = | sin 2 (7r<i/2). It confirms that type of heat conduction does not depend on 
concrete choice of on-site potential function. 

5 Conclusion 

We have considered the heat conduction process in the ID lattice of hard particles with periodic on-site 
potential. Analytical treatment predicts that for zero diameter of the particles the system will be completely 
integrable regardless the exact shape of the on-site potential. Therefore the heat conductivity will be infinite. 
For any nonzero size of the particles the heat transfer is governed by diffusion of quasiparticles giving rise to 
finite heat conductivity. The value of the heat conduction coefficient computed by the analytical treatment is 
in line with numerical simulation data. This coincidence is very profound if speaking about the asymptotic 
scaling behavior of the heat conduction coefficient in the cases of small and large particle sizes, as well as 
for the case of high temperatures. The characteristic behavior of the heat conduction coefficient does not 
depend on the exact shape of the on-site potential function. 

The above results mean that there exists a new class of universality of ID chain models with respect 
to their heat conductivity. The limit case of zero-size particles is integrable, but the slightest perturbation 
of this integrable case by introducing the nonzero size leads to drastic change of the behavior - it becomes 
diffusive and the heat conduction coefficient converges. It should be stressed that this class of universality, 
unlikely the systems with conserved momentum, cannot be revealed by sole numerical simulation. The 
reason is that the correlation length (as well as the heat conduction coefficient) diverges as the system 
approaches the integrable limit; therefore any finite capacity of the numerical installation will be exceeded. 
That is why the analytical approach is also necessary. 

The authors are grateful to Russian Foundation of Basic Research (grant 01-03-33122) and to RAS 
Commission for Support of Young Scientists (6th competition, grant No. 123) for financial support. O.V.G. 
is grateful to Russian Science Support Foundation for the financial support. 




(21) 



V(r) = lim V(e: r). 
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